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FOREWORD 

This  report  describes  an  algorithm  for  the  solution  of  the  incompressible  Navier-Stokes  equa¬ 
tions  in  which  the  convective  and  diffusive  terms  are  treated  independently.  This  treatment  allows 
for  the  use  of  a  second  order  Godunov  method  to  discretize  the  nonlinear  convection  terms.  The 
resulting  scheme  is  stable  for  high  Reynold’s  number  flows. 
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CHAPTER  1 
INTRODUCTION 

The  major  difficulty  in  designing  algorithms  for  problems  which  have  both  convective  and 
diffusive  terms  arises  from  the  fact  that  the  algorithm  should  work  well  for  both  parabolic  problems 
(when  diffusion  dominates)  as  well  as  for  hyperbolic  problems  (when  convection  dominates).  In  the 
diffusion  dominant  case,  standard  finite  element  or  finite  difference  methods  which  yield  space  cen¬ 
tered  discretizations,  coupled  with  a  stable  time  discretization,  will  generally  yield  accurate  results. 
However,  central  differencing  of  the  convection  terms  will  lead  to  an  unstable  (in  space)  scheme, 
evidenced  by  spurious  oscillations.  If  the  oscillations  are  local,  this  problem  can  be  alleviated  with 
an  appropriate  mesh  refinement.1  This  procedure  will  fail  if  either  the  oscillations  are  global,2  or 
they  appear  in  smooth  regions  of  the  flow,3  which  could  mislead  one  into  refining  the  grid  in  the 
wrong  area.  One  remedy  to  this  problem  is  to  incorporate  a  priori  information  about  the  solution 
into  the  difference  scheme,  which  results  in  schemes  which  are  nearly  upwind  when  convection 
dominates  and  are  nearly  centered  when  diffusion  dominates.4  Although  this  procedure  works  well 
for  one  dimensional  steady-state  problems,  its  generalization  to  multidimensional  problems  has  not 
yet  been  successful  due  to  the  complicated  behavior  of  the  solutions,5  and  the  excessive  diffusion 
resulting  when  the  operators  in  each  direction  are  treated  as  in  the  one  dimensional  case.6  An  alter¬ 
native  approach  is  to  split  the  processes  of  convection  and  diffusion,  and  use  independent  stable  and 
accurate  schemes  for  each  part. 

Various  studies  have  been  performed  concerning  "viscous  splitting"  in  References  7  through 
10.  When  both  the  diffusive  and  convective  steps  are  solved  exactly,  Beale  and  Majda7  proved  that 
the  error  is  O  (vt2)  in  L2,  where  x  =  A t  is  the  time  step,  and  v  is  the  viscosity  (diffusion 
coefficient).  Characteristic  Galerkin  methods,  which  may  be  considered  as  viscous  split  algorithms 
using  the  method  of  characteristics  for  the  convection  step,  applied  to  convection  diffusion  problems 
were  analyzed  in  References  8  and  9.  However,  both  treatments  were  only  first  order  accurate  in 
the  time  step,  due  to  the  use  of  the  backward  Euler  method  for  integrating  the  diffusion  step.  A 
second  order  viscous  split  method  was  analyzed  in  Reference  10  for  the  solution  to  steady  state  one 
dimensional  convection  diffusion  problems,  and  some  of  the  ideas  developed  in  that  paper  are  used 
here  as  well. 

This  report  is  concerned  with  applying  viscous  splitting  to  the  incompressible  Navier  Stokes 
equations,  employing  a  higher  order  Godunov  method  for  the  convection  step.  These  methods  have 
been  demonstrated  to  be  extremely  successful  for  solving  a  wide  variety  of  applications  requiring 
solutions  of  hyperbolic  systems  (see  e.g.,  References  11  and  12  and  the  references  cited  therein). 
We  remark  that  the  algorithm  of  Bell  et  al.,13  which  does  not  independently  split  the  diffusion  and 
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convection  terms,  is  similar  to  the  split  algorithm  employed  here  in  the  sense  that  both  employ  a 
Godunov  scheme  for  the  discretization  of  the  convective  terms  as  well  as  a  similar  projection 
scheme  for  the  treatment  of  the  incompressibility  constraint. 
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CHAPTER  2 
VISCOUS  SPLITTING 


The  time  dependent  incompressible  Navier-Stokes  equations  in  a  bounded  domain  fi  are 

u,  -  vAu  +  (u-V)u  =  -Vp  +  F,  (2-la) 


V  u  =  0, 


(2-lb) 


with  initial  conditions 


u(x,0)  =  u0(x),  for  x  e  Q, 

and  boundary  conditions 

B(u)  =  b(x)  for  x  e  dQ, 


(2-lc) 


(2- Id) 


where  u  =  (u  ,v )  is  the  velocity,  p  is  the  pressure,  v  is  the  viscosity,  and  F  is  a  source  function. 
Equation  (2-lb)  is  the  incompressibility  constraint  and  (2- Id)  specifies  boundary  conditions  of  either 
Dirichlet  and/or  Neumann  type  on  various  pans  of  the  boundary  9Q.  When  the  problem  is  non- 
dimensionalized  by  a  characteristic  length  and  velocity  the  Reynolds  number  Re  is  defined  by 


v 


The  solution  of  (2-1)  will  be  approximated  by  solving  the  sub-problems  of  diffusion  and  con¬ 
vection  separately.  That  is,  (2-1)  will  be  split  into  Stokes  equation  (diffusion) 

w,  -  vAw  =  -Vp  +  F, 

V  w  =  0, 


w(x,0)  =  w0(x),  for  a  e  Q, 

B(w)  =  b,  for  x  g  9f2,  (2-2) 

and  a  convection  step 

z,  +  (Z'V)z  =  0, 

z(x,0)  =  Zq(x),  for  X  G  Q,  (2-3) 


where  the  treatment  of  the  boundary  conditions  for  the  convection  step  will  be  discussed  later.  Let 
S(t)  denote  the  solution  operator  to  (2-2)  after  a  timestep  x.  That  is,  w  =  S(z)w0  is  the  solution  to 
(2-2)  at  time  t  =  x.  Also,  let  C(x)  be  such  that  z  =  C(x)z0  is  a  solution  to  (2-3)  at  time  t  =  x  .  The 
viscous  split  time  stepping  procedure  for  the  approximate  solution  to  (2-1)  is 
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un  =  S(t/2)C<t)S(t/2)u„_i 


(2-4) 


where  un  is  an  approximation  to  u(x,«x).  This  type  of  splitting  is  often  referred  to  as  Strang  split¬ 
ting.14  In  addition,  consider  a  split  algorithm  in  which  the  diffusion  operator  is  replaced  by  its 
Crank-Nicolson  approximation  S  (trapezoid  rule  in  time).  Then  w  =  S(t)w0  solves  the  problem 


X  X  *  * 

-VyAw  +  W  =  VyAw0  +  w0  -  xVp  +  tF  , 


V  w  =  0,  for  x  €  Q, 

B(w)  =  b(x),  for  x  e  BQ,  (2-5) 


where  the  superscript 
becomes 


indicates  an  approximation  at  time  — . 


In  this  case  the  splitting  algorithm 


u„  =  S(t/2)C(t)S(t/2)u„_, 


1" 


[S(x/2)C(T)S(t/2)J  uc 


(2-6) 


2-2 


NSWC  TR  89-136 


CHAPTER  3 

SPATIAL  DISCRETIZATION 

As  mentioned  in  the  introduction,  the  main  advantage  of  viscous  splitting  is  that  different  and 
independent  discretizations  of  the  diffusion  and  convection  steps  may  be  used.  In  particular,  it 
allows  for  the  use  of  schemes  of  the  type  discussed  in  References  11  and  12,  which  have  proven  to 
be  very  successful  for  the  solution  of  hyperbolic  conservation  laws. 

In  this  report,  a  second  order  Godunov  method  with  an  approximate  Riemann  solver  is 
employed  for  the  discretization  of  the  convective  step.  For  the  two  dimensional  hyperbolic  system 

u/  +  F*  +  Gy  =  0  . 

This  algorithm  is  outlined  in  the  following  steps: 

1.  Compute  and  limit  the  slopes  in  the  x  and  y  directions  for  each  component  of  the  system.  This 
produces  a  piecewise  linear  (not  necessarily  continuous)  profile  for  u  in  each  direction.  The 
slope  limiting  is  performed  so  that  this  piecewise  linear  approximation  does  not  introduce 
spurious  extrema. 

2.  Compute  predicted  values  at  the  half  time  step.  The  predicted  values  are  computed  by 
differencing  F  and  G  at  cell  centers  using  values  for  u  at  cell  edges  from  the  linear  profiles 
generated  from  the  previous  step.  The  predicted  values  are  made  approximately  divergence 
free  by  adding  the  gradient  of  the  pressure  computed  from  the  previous  Stokes  solve.  This 
modification  is  necessary  to  preserve  the  formal  second  order  accuracy  of  the  viscous  split  rou¬ 
tine. 

3.  Compute  numerical  fluxes  along  each  edge  by  approximately  solving  a  Riemann  problem  in 
each  direction  based  on  values  from  the  prediction  step.  These  fluxes  must  be  consistent  with 
the  physical  flux,  i.e., 

F*  (u,u)  =  F(u), 

where  F/,(u,v)  is  the  numerical  flux  with  states  u  and  v  at  the  left  and  right,  respectively. 

4.  Correct  the  solution  using  the  computed  numerical  fluxes. 

The  particular  implementation  employed  for  the  results  in  this  report  is  described  in  detail  by 
Davis.11  The  slope  limiting  performed  with  these  schemes  can  be  shown  to  have  the  same  effect  as 
adding  artificial  viscosity  only  in  regions  where  the  concavity  of  the  solution  is  charging  abruptly, 
such  as  near  shocks  or  interior  or  boundary  layers.  The  solution  of  the  Riemann  problem  associated 
with  the  convection  terms  of  (2-3)  results  in  an  upwind  scheme  depending  on  the  sign  of  the 
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velocity  component.  Along  Dirichlet  boundaries  only  the  normal  component  of  the  velocity  is 
prescribed,  which  is  the  usual  condition  for  hyperbolic  problems.  The  Courant-Friedrichs-Levy  sta¬ 
bility  condition  restricts  the  timestep  by  cxlh  =  X  <  1  ,  where  c  =  max(  I  2 ui}  I  +  I  2vtJ  I  ).  This 

restriction  can  be  improved  by  a  factor  of  two  (c  reduced  by  one  half)  using  a  non-conservative 
scheme,  and  with  further  modifications  can  be  based  on  the  maximum  component  instead  of  the 
sum  (see  e.g.,  Reference  13). 

The  Stokes  solver  of  the  Galerkin  Finite  Difference  Method  (GFDM)  3,5,15,16  was  used  for  the 
discretization  of  S.  We  outline  this  method  for  the  case  of  Dirichlet  data.  Further  details  as  well  as 
details  of  the  appropriate  modifications  for  more  general  boundary  conditions  may  be  found  in 
Reference  3. 


This  method  begins  with  an  approximation  of  divergence  at  cell  centers  given  by 

U)(  +  l/2,;+l/2  —  2h  (W,+h)  +  l  _  Ut.j+l  U‘  +  Lj  ~  ui.j  +  V.+l,/  +  l  ~  Vi+l,j  +  Vi.j  +  )  ~  Vi.j  j 
and  an  approximation  to  the  gradient  defined  at  cell  vertices  by 


(G*< »ij 


2  h 


|^i+l/2,;+l/2  "  0i-l/2.)+l/2  +  4>i+V2.j-V2  ~  $1- 1/2,;' -1/2 
'Qi+m,j+V2  -  Qi+VZj-m  +  $i-V2,j+\a  ~  ^i-l/2J-V2 


(3-2) 


for  a  scalar  function  0  defined  at  cell  centers.  The  linear  operator  maps  a  space  of  mesh  vectors 
defined  at  the  vertices  to  a  space  of  mesh  scalars  Wh  defined  at  cell  centers,  and  Gh  maps  Wh 
into  \h.  In  Reference  15  these  spaces  are  equipped  with  standard  l2  inner  products  such  that 
D*  =  Ga*  holds,  where  G^  denotes  the  adjoint  of  Gh .  We  remark  that  these  definitions  for  Dh  and 
Gh  are  the  same  discretizations  that  arise  from  a  finite  element  approximation  with  bilinear  veloci¬ 
ties  and  piecewise  constant  pressures.  An  important  observation  employed  by  the  GFDM  is  that 
Dh  =  {'j'*+1/2’/+1/2:  }  form  a  basis  for  the  kemal  of  where 


i+l/2,;+l/2 


,  ((-l)/--/,(-l)*"'+1/  for  i-k  ,k+\;j=l  ,1+1 
(0,0)  for  all  other  i  J 


(3-3) 


The  GFDM  approximation  to  (2-5)  is  then  to  find  wh  =  £a(y¥J+1/2j+1/2  such  that 


-v-jaaw>.  +  "*.¥ 


v-jA*wo  ,h  +  w0  ,h  +  xF*,¥ 


for  each  ¥  e  Dh  , 


(3-4) 


where  Aa  denotes  for  example  the  standard  5  point  discrete  Laplacian,  and  |vj  denotes  an  inner 
product  on  \h .  The  values  for  <xl}  along  Dirichlet  boundaries  are  determined  by  the  prescribed  data 
provided  certain  compatibility  conditions  are  met  (again  see  References  3  or  15).  The  pressure  is 
eliminated  from  this  computation  due  to  the  fact  that  j-GA/?*,¥j  =  jp*,DA¥j  =0  for  each 
¥  e  Dh.  However,  Ghp*  can  be  computed,  for  use  in  the  convection  (predictor)  step,  using  the 
discretization  of  (2-5),  after  wh  is  determined  from  (3-4). 
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CHAPTER  4 
NUMERICAL  RESULTS 

Since  the  argument  for  order  of  accuracy  was  only  formal,  we  first  conduct  a  numerical  con¬ 
vergence  study  on  a  smooth  time  dependent  problem.  The  domain  for  this  problem  is  the  unit 
square  with  homogeneous  Dirichlet  data  and  initial  conditions 

u  ( x  ,y  ,0)  =  0.5sin(7ty  )cos(ny  )sin2(Tr ) 

v  (x  ,y  ,0)  =  -0.5sin(7Lt  )cos(ju  )sin2(7ty )  . 

Solutions  were  computed  on  a  sequence  of  uniform  grids  with  N  by  N  cells  with  N  =  4,  8,  ....  64. 
The  solutions  were  computed  to  time  T  =  0.125,  with  the  timestep  i  =  h/2  which  corresponds  to  the 
Courant  number  X  =  .323.  Listed  in  Table  4-1  are  the  /2  differences  EN  =  ||uA  -  ||/2  ,  where 

the  1 2  norm  is  computed  on  the  coarser  N  by  N  grid.  From  this  table  the  rates  of  convergence  for 
Re=0  and  Re=1000  are  observed  to  be  at  least  two.  For  Re=  the  rate  is  slightly  less  although  the 
asymptotic  rate  has  not  been  attained  for  this  case. 

TABLE  4-1.  1 2  DIFFERENCES  AND  CONVERGENCE  RATES 


Re=0  (Stokes) 

Re=1000 

Re=°o  (Euler) 

N-2N 

En  Rate 

£\  Rate 

£,v  Rate 

4-8 

,848e-2 

,755e-2 

,798e-2 

2.22 

1.46 

1.45 

8-16 

. 183e-2 

,275e-2 

.29 1  e-2 

2.09 

2.36 

2.26 

16-32 

,429e-3 

,534e-3 

,606e-3 

2.03 

2.02 

1.59 

32-64 

. 105e-3 

. 1 32e-3 

,201e-3 

The  algorithm  is  further  demonstrated  on  a  problem  of  2-D  channel  flow  over  a  full  step  at 
Re=100  (based  on  the  inlet  height  and  average  inlet  velocity).  For  this  problem,  the  inlet  velocity  is 
parabolic,  the  channel  has  height  2  and  length  30,  the  step  is  1  unit  high  and  2  units  wide,  and  the 
front  face  of  the  step  is  6  units  from  the  inlet.  Figure  4-1  displays  the  streamline  contours  of  the 
steady  state  solutions  obtained  using  the  viscous  split  algorithm  and  the  unsplit  GFDM16  using 
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conservative  centered  differencing  for  the  convective  terms,  on  identical  30  by  16  grids.  (The  plots 
are  stretched  vertically  by  a  factor  of  two  for  clarity.)  In  this  case  the  viscous  split  solution  produces 
results  similar  to  streamline  upwind  type  schemes,  as  described  in  References  15  through  17,  in  the 
sense  that  spurious  oscillations  near  the  front  face  of  the  step  are  eliminated  and  the  length  of  the 
recirculation  region  behind  the  step  is  not  decreased  by  excessive  artificial  diffusion. 


VISCOUS  SPLIT  GODUNOV  SCHEME 


FIGURE  4-1.  COMPUTED  STREAMFUNCTION  CONTOURS  OF  STEADY  STATE 
SOLUTIONS  AT  P.e=100  ON  A  30  BY  16  GRID 


Further  tests  on  steady  state  problems  for  values  of  Re  <  1000  were  performed  for  flow  over  a 
rearward  facing  step  and  compared  to  the  experimental  results  of  of  Reference  18  (as  was  done  for 
the  unsplit  GFDM  in  Reference  16)  with  great  success,  but  these  results  are  not  shown  here. 
Instead,  the  time  evolution  of  a  flowfield  at  Re= 10000  with  initial  velocities  determined  from  the 
steady  state  Stokes  solution  is  displayed  in  Figure  4-2.  Although  the  120  by  32  grid  used  in  the 
calculation  is  insufficient  to  resolve  all  details  of  the  flow,  particularly  near  the  step,  we  believe  a 
reasonable  depiction  of  the  gross  properties  of  the  flow  are  preserved.  The  multiple  recirculation 
regions  have  also  been  observed  experimentally,18  as  well  as  for  steady  state  computations  at  lower 
Reynolds  numbers.16 

The  generation  of  stable  (but  not  overdiffused)  approximations  on  grids  where  all  features  of 
the  flow  are  not  resolved  has  important  implications  for  developing  an  adaptive  mesh  refinement 
strategy,  as  demonstrated  in  References  2  and  10  for  one  dimensional  steady  problems.  This 
feature,  in  addition  to  the  property  that  an  implicit  time  integration  of  the  parabolic  terms  results  in 
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a  positive-definite  symmetric  system,  makes  viscous  splitting  an  attractive  approach  for  solving 
Navier-Stokes  as  well  as  other  classes  of  problems  with  convective  and  diffusive  terms. 


T=39.49 


FIGURE  4-2.  TIME  EVOLUTION  OF  COMPUTED  STREAMFUNCTION  CONTOURS 
AT  Re=  10000  USING  VISCOUS  SPLITTING  ON  A  120  BY  32  GRID 
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